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Abstract 

The  purpose  of  this  study  was  to  use  the  computer 
technique  of  atomistic  simulation  to  investigate  the  equation 
of  state  for  copper.  Specifically,  this  study  used  a 
deterministic  technique  called  molecular  dynamics  with  a 
relatively  new  method  to  determine  interatomic  forces  called 
the  embedded-atom  method.  The  basic  approach  was  to  modify  a 
computer  program  which  successfully  works  on  the  Cray 
supercomputer  to  work  on  other  computers  which  are  more 
readily  available.  As  a  final  step,  the  program  was  modified 
to  work  on  a  parallel  processing  computer  system  in  order  to 
study  the  time  savings  made  possible  by  solving  parts  of  the 
problem  in  parallel. 

This  study  found  that  the  molecular  dynamics  approach  and 
the  embedded-atom  method  gave  an  accurate  prediction  of  the 
copper  lattice  constant  at  room  temperature  and  of  thermal 
expansion  at  higher  temperatures.  In  order  to  determine  the 
melting  temperature,  the  simulated  crystal  was  then  heated 
past  its  melt  point,  cooled  until  it  froze,  and  reheated  until 
it  melted  again.  The  experiments  found  a  transition 
temperature  very  near  1350°  K,  which  is  close  to  the  accepted 
melt  point  of  1357.6°  K. 


The  experiments  with  the  parallel  processing  system  showed 
that  a  program  running  in  a  parallel  mode  could  realize  time 
savings  of  a  factor  of  3.75  over  serial  processing  systems. 
Because  the  code  was  essentially  modular  in  nature,  the  most 
extensive  modifications  to  make  it  run  in  parallel  involved 
separating  the  necessary  subroutines  and  adding  the 
communications  routines  so  data  could  be  passed  back  and  forth 
during  the  calculations.  The  modifications  to  run  in  parallel 
did  not  make  the  program  more  difficult  to  manipulate  and 
could  quickly  be  removed  in  order  to  run  the  program  on  a 
normal,  serial  computer. 


A  STUDY  OF  THE  EQUATION  OF  STATE  FOR  COPPER 


I .  INTRODUCTION 

Atomistic  simulation  is  a  computer  technique  used  to  study 
matter.  The  general  idea  is  for  a  computer  to  simulate  a 
system  of  discrete  atoms.  As  such,  atomistic  simulation 
differs  from  continum  methods  such  as  fluid  dynamics.  The 
atoms  are  allowed  to  react  to  random  motions  and  to  interatomic 
forces,  and  the  computer  determines  each  atom's  position  and 
velocity  at  discrete  time  intervals.  The  computer  then 
calculates  properties  such  as  temperature  and  pressure,  which 
are  functions  of  atomic  position  and  velocity.  These  functions 
are  averaged  over  several  thousand  time  steps  so  that  the 
overall  state  of  the  atomic  system  can  be  determined. 

An  interatomic  potential,  the  energy  between  atoms  as  a 
function  of  distance,  is  a  way  to  express  how  atoms  interact 
with  each  other.  The  forces  on  each  atom  can  be  calculated 
after  the  total  energy  of  a  system  has  been  calculated  (15:1). 
Atomistic  simulation  has  the  advantage  of  being  able  to  work 
with  many  different  interatomic  potential  functions.  On  the 
other  hand,  the  application  of  the  results  of  atomistic 
simulation  experiments  depend  heavily  upon  the  accuracy  of  the 
interatomic  potential  function  applied  within  the  computer 


program. 


Atomistic  simulation  gives  the  experimenter  complete 
control  of  the  external  environment  which,  in  turn,  allows 
environments  that  may  be  unavailable  in  the  laboratory.  The 
computer  can  generate  large  quantities  of  data  on  the  material 
studied,  and  extremely  short-lived  phenomena  can  be  simulated, 
such  as  the  motion  of  phonons  within  a  crystal  (4)  and  the 
melting  of  metallic  crystals  (8). 

The  purpose  of  this  study  was  twofold.  First,  the  equation 
of  state  for  copper  was  investigated  by  using:  1)  a  particular 
atomistic  simulation  technique,  molecular  dynamics,  and  2)  a 
relatively  new  technique  to  determine  interatomic  force,  the 
embedded  atom  method.  Second,  this  study  used  both  serial  and 
parallel  processing  computers  in  an  attempt  to  generate 
accurate  data  in  the  least  amount  of  time  and  to  validate  the 
use  of  parallel  processing  techniques  on  this  particular  type 
of  problem. 

General  Approach 

This  series  of  experiments  used  a  computer  code  developed 
by  Baskes  and  Foiles  (1)  for  use  on  the  Cray  computer  and 
modified  by  the  researcher  to  run  on  the  ELIXI  System  6400 
under  the  EMBOS  11.3  operating  system.  The  modified  code  was 
run  extensively  on  the  ELIXI  system  in  order  to  validate  the 
output  and  to  perform  several  computer  experiments  on  both 
crystallized  and  molten  copper. 


The  code  was  further  modified  to  run  on  the  Intel  Personal 
Supercomputer  (iPSC).  The  iPSC  is  a  computer  which  employs  up 


to  thirty-two  concurrent  processors.  Each  processor  can  work 
independently  of  the  others.  Ideally,  any  given  program  is 
solved  much  more  quickly  than  on  a  conventional  serial 
processor.  Unfortunately,  not  every  aspect  of  the  molecular 
dynamics  program  lends  itself  to  a  parallel  processing 
approach.  The  techniques  and  results  are  discussed  later  in 
this  text. 

The  steps  which  the  computer  uses  to  solve  the  problem  are: 
1)  establish  matrices  that  describe  the  type  of  material, 
position,  and  velocity  of  every  atom  in  the  crystal  to  be 
studied  in  the  experiment,  2)  determine  the  energies  and  forces 
on  each  atom  for  the  current  time  step,  3)  use  those  forces  to 
project  atomic  accelerations  and  displacements  over  the  time 
step,  4)  use  the  displacements  and  forces  to  provide  snapshots 
of  the  locations  of  the  atoms  at  given  time  intervals,  5) 
calculate  state  functions  such  as  system  energy,  internal 
pressure,  internal  temperature,  and  stress,  and  6)  begin  the 
next  time  step  by  calculating  the  forces  and  continuing  with 
the  next  displacements.  When  a  specified  number  of  time  steps 
is  complete,  averages  taken  over  all  time  steps  (after  a  short 
period  to  let  forces  even  out)  are  used  to  determine  the  final 
state  of  the  original  material.  The  final  state  can  then  be 
stored  in  a  computer  file  and  used  as  initial  positions  for 
subsequent  computer  experiments. 


Significance  of  this  Study 

This  study  is  significant  in  three  areas:  1 )  as  a 
straightforward  study  of  solid  and  liquid  copper,  2 )  as  a 
validation  of  a  particular  embedded-atom  potential  for  copper, 
and  3)  as  an  expansion  in  the  use  of  parallel  processing 
computers . 

Copper  is  well  known  for  high  electrical  conductivity  and 


is 


for  efficient  heat  diffusion.  Any  further  study  of  metals  in 
general  and  of  copper  in  particular  can  add  further  insight 
into  the  nature  of  matter.  As  new  alloying  techniques  and  new 
methods  to  alter  surface  molecular  arrangements  continue  to 
develop,  extensive  knowledge  of  material  will  become 
increasingly  critical  (16:227,228). 

Should  molecular  dynamics  become  a  more  necessary  scien¬ 
tific  tool,  programs  must  be  applied  to  computers  which  are 
more  readily  available  than  a  Cray  supercomputer.  During  this 
study  the  modified  code  operated  on  other  computers  which  are 
less  costly  and  more  accessible.  This  advance  indicates  that 
molecular  dynamics  can  be  used  by  students  and  other 
investigators  more  easily  than  before. 

This  study  also  shows  that  the  embedded-atom  method  gives 
useful  interatomic  potentials.  The  other  advantages  of  the 
embedded  atom  method,  such  as  its  applicability  to  other  metals 
and  alloys,  will  make  this  approach  useful  in  other  atomistic 
simulation  approaches. 

Finally,  the  Air  Force  is  constantly  attempting  to  increase 
computer  control  of  aircraft  and  other  weapons  systems. 
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Parallel  processing  promises  to  be  much  faster  than 
conventional  serial  processing,  and  that  speed  will  be 
extremely  useful  in  managing  interdependent  tasks  often  found 
in  flight  control,  communications,  and  electronic  warfare 
systems.  Advances  in  the  application  and  understanding  of 
parallel  processing  are  important  for  both  scientific 
investigation  and  for  future  Air  Force  needs. 

Order  of  Presentation 

The  specifics  of  molecular  dynamics  are  detailed  in  Chapter 
Two.  The  discussion  includes  the  important  calculations 
completed  during  the  simulation  and  the  different  methods  used 
to  determine  the  interatomic  potentials  needed  for  force 
calculations. 

Chapter  Three  discusses  parallel  processing,  how  it  differs 
from  conventional  serial  processing,  and  how  molecular  dynamics 
and  the  embedded  atom  method  can  be  adapted  to  parallel 
processing  techniques. 

Chapter  Four  discusses  the  computer  experiments  and  their 
results.  The  ELIXI,  which  uses  conventional  serial  computing, 
gave  extensive  data  on  the  equation  of  state  of  copper.  The 
iPSC  parallel  processing  computer  was  used  to  show  the 
potential  time  savings  by  using  parallel  computing  techniques. 

The  last  chapter.  Chapter  Five,  discusses  conclusions  and 
makes  recommendations. 


II.  ATOMISTIC  SIMULATION 


The  advent  of  the  programmable  electronic  computer  has 
made  atomistic  simulation  a  feasible  experimental  tool. 
Physicists  in  the  laboratory  can  only  observe  statistical 
parameters  such  as  temperature  and  pressure.  Although  the  sum 
of  the  forces  on  each  atom  determines  those  observable 
parameters ,  no  one  could  observe  the  motions  of  each 
individual  atom.  Atomistic  simulation  makes  the  forces  on 
each  atom  observable  and  the  computer  can  use  the  sum  of  the 
forces  on  each  atom  in  order  to  calculate  the  parameters 
involved  in  the  equation  of  state. 

The  simulation  technique  used  in  this  study  is  called 
molecular  dynamics.  Molecular  dynamics  is  almost  completely 
deterministic.  After  random  numbers  with  a  Boltzmann 
distribution  are  used  to  assign  initial  velocities  to  each 
atom,  atomic  motions  depend  entirely  on  interatomic  potential 
calculations.  The  computer  calculates  subsequent  atomic 
motions  and  uses  those  motions  to  determine  parameters  of 
state. 

The  purpose  of  this  chapter  is  to  discuss  the  mechanics  of 
atomistic  simulation.  The  opening  section  discusses  how  the 
lattice  is  simulated  within  the  computer  and  following 
sections  discuss  subsequent  steps  of  the  solution  process: 
calculating  the  interatomic  forces,  integrating  those  forces 
and  determining  the  atomic  motions,  and  calculating  the 
parameters  of  state  from  the  atomic  motions.  This  chapter 


concludes  by  discussing  other  pertinent  problems  such  as  the 
determination  of  the  length  and  number  of  time  steps  necessary 
to  solve  a  given  problem,  what  alternative  configurations  can 
be  computed  from  the  molecular  dynamics  approach,  and  what 
techniques  are  available  to  save  computer  time. 

Building  the  Lattice 

A  crystal  of  256  copper  atoms  is  used  for  all  computer 
runs  in  this  study  because  it  provides  a  good  compromise 
between  the  desire  for  a  large  crystal  and  the  limitations  on 
computer  speed  and  capacity  (13:20).  Copper  forms 
face-centered  cubic  (fee)  crystals  and  each  unit  cell  of  an 
fee  system  contains  four  atoms  < 19 : 12 ; 12 : 6 ) .  Figure  2-1  shows 
how  the  atoms  are  arranged  in  a  single  fee  crystal.  A  crystal 
of  four  lattice  constants  on  each  side  gives  (4  x  4  x  4  =  64) 
unit  cells  or  (64  x  4  =  256)  atoms  in  all. 

Each  atom  in  the  crystal  is  then  described  by  its  position 
and  its  velocity  in  each  of  the  three  spatial  dimensions. 

Hence,  each  atom  must  have  six  coordinates: 

spatial  position  (x,  y,  z)  and  velocity  (dx  dy  dz ) 

(dt,  dt,  dt) 

Initial  positions  are  determined  by  the  type  of  lattice 
structure  and  by  the  limits  on  crystal  size.  Velocities  are 
generated  by  assigning  random  numbers  with  a  Boltzmann 
distribution  based  on  the  initial  temperature. 

After  the  matrices  to  describe  the  crystal  are  built  and 
used  in  a  computer  experiment,  their  positions  and  velocities 
can  be  stored  in  a  file.  Subsequent  computer  runs  can  use  that 

2-2 


M 


<“?V 


m 


.  v  \rv 


Figure  2-1.  A  Face-Centered  Cubic  Crystal 

(from  12:4) 
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file  as  initial  positions  and  velocities  so  the  lattice  does 
not  need  to  be  rebuilt.  With  this  procedure,  for  example,  a 
single  crystal  can  be  heated  and  cooled  over  the  course  of 
several  computer  runs. 


Boundary  Effects 

In  order  to  eliminate  boundary  effects,  the  atomic 
system  is  simulated  to  be  surrounded  by  images  of  itself. 

This  is  illustrated  on  Figure  2-2.  The  distances  between  the 
actual  crystal  and  its  images,  the  periodic  bounds,  are 
specified  by  the  program  operator.  If  an  atom  crosses  a 
periodic  boundary  during  any  time  step,  that  atom  is  placed  on 
the  other  side  of  the  crystal  as  shown  in  Figure  2-3. 
Experience  with  the  program  demonstrated  that  if  the  operator 
imposes  unrealistic  periodic  bounds,  the  simulation  responds 
with  unrealistic  pressures  and  temperatures. 

An  alternative  technique  to  minimize  boundary  effects  on  a 
single  crystal  would  be  to  use  a  very  large  crystal  with 
vacuum  boundaries.  Surface  area,  in  general,  tends  to  expand 
at  a  rate  proportional  to  the  square  of  any  increase  in  the 
number  of  atoms;  mass  and  volume  increase  in  proportion  to  the 
cube  of  any  increase.  As  a  crystal  becomes  larger,  therefore, 
the  proportion  of  surface  area  to  mass  would  become  smaller, 
and  surface  effects  would  decrease.  A  crystal  of  one  million 
atoms  or  more,  for  example,  would  have  such  minimal  boundary 
effects.  Extremely  large  crystals  are  not  practical  for  the 
computers  available  at  this  time. 


Periodic  Boundaries 
B  interacts  with  the  image  of  A 
(arrows  indicate  velocity  vectors) 

(from  17:11) 
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Force  Calculations 

Force  and  energy  calculations  are  the  core  of  any 
molecular  dynamics  program.  The  general  idea  is  that  the 
forces  on  an  atomic  system  can  be  computed  by  summing  the 
forces  on  each  individual  atom: 


U  (1,  N)  = 


DUilrJ  +  £  n2(r±,  r3 )  +  £  u3  (rif  rj(  rk) 

i  i,j  i,j,k 

i<j  i<j<k 


+  ...  +  UN  (  Ti  ,  r  2  ,  1 3  t  ...  £"r-k  ) 


(2.1) 


where  Ui(r±)  is  the  force  on  atom  i  imposed  by  the  external 
environment,  u2(ri.,  r., )  is  the  intermolecular  force  between 
the  two  atoms  i  and  j ,  u3  is  the  three-body  force  between  i , 
j,  and  k,  and  so  on  up  to  uN  which  considers  the  entire  system 
of  N  atoms  interacting  together.  This  approach  obviously 
becomes  unwieldy  for  a  large  system,  so  computer  simulations 
must  use  less  complex  approximations  (6:51). 

The  advantage  of  the  pair  potential  method  is  that  it  is 
simple  and  straightforward  (13:7).  This  method  assumes  that 
all  forces  beyond  the  second  order  terms  are  negligible.  The 
overall  force  on  each  atom  is,  therefore,  computed  as  a 
superposition  of  the  forces  on  each  two  atom  pair  (6).  The 
computer  uses  a  truncation  of  equation  (2.1),  where,  with  no 
external  forces: 


U  (1,  N)  =  lu2(rl(  r3) 


(2.2) 
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The  pair  potential  method  necessarily  assumes  spherical 
symmetry,  and  has  proven  successful  in  simulations  of  noble 
gasses  (13:7).  Where  electrons  are  less  tightly  bound,  such 
as  in  metals,  however,  the  pair  potential  method  does  not  give 
satisfactory  results  and  the  force  calculations  must  become 
more  complicated. 


The  Embedded-Atom  Method  (EAM) 


Because  electrons  in  metals  are  less  tightly  bound 
than  electrons  in  the  noble  gasses,  the  atomic  environment  can 
be  modeled  as  nuclear  cores  surrounded  by  both  their  own 
electrons  and  by  unbound  electrons  acting  like  an  electron 
gas.  The  basic  approach  of  the  embedded-atom  method  is  to 
calculate  the  force  required  to  embed  each  atom  in  a  host 
lattice  of  all  the  other  atoms  in  the  crystal  (2:6443).  Local 
electron  density  is  computed  from  the  electron  density  of  the 
atom  and  the  background  density  of  the  electron  gas  in  the 
lattice  (8:7984).  Pair  potentials  are  used  to  determine  the 
repulsion  between  atomic  cores ,  and  the  total  force  on  each 
atom  is  the  sum  of  the  core-core  repulsion  and  the  force 
imposed  by  the  unbound  electrons  (4:697). 

The  general  embedded-atom  method  is  not  significantly  more 
complicated  to  use  than  pair  potentials,  and  the  actual 
calculations  are  not  especially  more  computer  intensive.  In 
addition,  the  embedded-atom  method  is  successful  at  describing 


more  complex  situations  and  is  especially  useful  for  studies 
of  alloys  as  well  as  for  pure  metals  (8:7984). 

Daw  and  Baskes  first  suggested  the  embedded-atom  method 
(2),  and  it  has  been  successfully  applied  to  hydrogen 
embrittlement  in  metals  (3),  phonons  in  transition  metals  (4), 
liquid  transition  metals  (7),  and  to  defects  in  metals  (2). 

The  success  of  the  embedded-atom  method  is,  in  the  final 
analysis,  what  justifies  its  use. 

The  first  task  of  an  embedded-atom  algorithm  is  to 
determine  the  electron  density  at  the  location  of  each  atom. 
This  density  is  approximated  by  a  superposition  of  atomic 
densities : 

P*,.  ±  =  L  PS  (Ri.a)  (2.3) 

j  (*i) 

where  /)h<1  is  the  host  electron  density  at  atom  i  imposed  by 
the  other  atoms  in  the  lattice,  A,“(R)  is  the  electron  density 
contributed  by  atom  j,  and  R±:1  is  the  distance  between  atoms 
i  and  j  (7:3410)  . 

A  known  electron  density  can  be  applied  to  an  embedding 
function  Fe±(p)  which  gives  the  energy  needed  to  embed  atom  i 
into  the  background  density  .  A  second  function,  a.a(Ri.j), 
is  simply  a  pair  potential  which  gives  the  core-core  pair 
repulsion  between  atoms  i  and  j . 

The  total  system  energy  is  then  given  by: 

1 

Etot  =  JCFSj.  (ph,  i)  t  2  2  (Ri  j  )  (2.4) 

2  i  j(*i) 


(7:3409) 


Values  for  the  embedding  function  Fe  and  the  core -core 
repulsion  function  <f>  have  been  obtained  empirically  by  fitting 


physical  properties.  Daw  and  Baskes  determined  Fe  and  by 
using  lattice  constants,  elastic  constants,  vacancy-formation 
energies,  and  sublimation  energies  in  fee  and  bcc  phases 
(7:3410).  In  the  computer,  the  values  are  stored  in  a  table 
and  splines  are  used  to  interpolate  the  specific  values  needed 
for  the  calculations.  The  file  which  contained  the  values  for 
the  copper  interaction  functions  was  provided  with  the 
original  computer  code. 

Integration  and  Computation 

Because  each  atom  exerts  force  on  its  neighboring  atoms, 
Newton's  Second  Equation  of  Motion: 

ma  =  F  (2.5) 

is  rewritten: 


d  d 

mi - (vA)  =  F±  =  -  [u  (i)]  (2.6) 

d  t  d  ri;) 


where  Fx  is  the  force  on  atom  i,  a  is  acceleration,  m  is  the 
mass  of  each  atom,  r13  is  the  distance  between  atoms  i  and  j, 
and  u(i)  is  the  interatomic  potential  on  atom  i.  Forces 
caused  by  the  interatomic  potential  determine  each  atom's 
acceleration  and  displacement,  the  equations  of  motion  are 
solved  over  time  by  using  a  fifth-order  predictor-corrector 
algorithm.  The  forces  are  then  summed  and  observable 
parameters  are  calculated. 


Calculation  of  Parameters 

All  state  functions,  such  as  temperature  and  potential 
energy,  are  given  by  algorithms  well  established  in  molecular 
dynamics  (10: 10; 13: 28, 29; 17: 10) . 

The  stress  tensor  is  of  particular  importance  because  it 
determines  internal  pressure,  and  that  internal  pressure 
influences  atomic  motions  in  subsequent  time  steps.  The 
stress  tensor,  which  describes  the  stress  within  the  crystal, 
is  determined  by  the  motions  of  each  atom  multiplied  by  its 
weight  and  reduced  by  the  interatomic  attractions  between  the 
atoms.  The  result  is  a  three-by- three  stress  tensor  which  is 
symmetric  to  the  diagonal  (1).  Readers  especially  interested 
in  the  stress  calculation  should  refer  directly  to  the  Baskes 
and  Foiles  computer  program. 

The  potential  energy  of  each  atom  is  determined  during  the 
force  calculations.  System  potential  energy  is  a  summation  of 
the  potential  energy  of  each  atom: 

U  =  Iu(i)  (2.7) 

i=l 

Kinetic  energy  is  calculated  by  a  straightforward: 

1  N 

K.E.  =  ---  m  L  v2  (2.8) 

2  i=l 

summed  over  each  atom's  calculated  velocity. 

Temperature  is  then  calculated  by: 

const  x  K.E. 

Temperature  =  -  (2.9) 

number  of  atoms 


where  the  constant  converts  from  units  of  energy  to  degrees 
Kelvin. 


WWfWWWWV»mfWUllHWUI 


i 


I 


► 


1  M 

<  P  >  =  —  Z  P  (tm)  (2.13) 

M  m=l 

Determining  the  Time  Steps  and  Length  of  Runs 

The  optimum  time  step  for  each  integration  is  a  compromise 
between  short  steps  for  accuracy  and  long  steps  which  save 
computer  time  and  still  generate  meaningful  data.  A  time  step 
of  0.001  x  10-12  second  (0.001  picosecond  or  10-15  sec)  gives 
the  required  compromise  (17:16). 

During  each  run  to  generate  data  on  the  reactions  of  a 
copper  crystal,  the  results  from  the  first  300  time  steps  (0.3 
ps)  are  not  added  into  the  final  averages  in  order  to  let  the 
crystal  react  to  the  changes  in  stress  imposed  by  the 
experiment.  This  delay  ensures  that  the  entire  crystal  is 
equilibrated  with  a  Maxwell-Boltzmann  distribution  of 
velocities  based  on  temperature.  The  crystal  is  then  left 
undisturbed  for  two  or  five  ps,  depending  on  the  nature  of  the 
particular  computer  run.  Both  instantaneous  and  average 
values  of  important  parameters,  such  as  temperature  and 
pressure,  are  printed  at  intervals  selected  by  the  operator  in 
order  to  show  trends  and  to  establish  that  there  are  no 
discontinuous  jumps  in  the  data  which  would  indicate 
inaccuracies  in  calculations  or  the  use  of  excessively  long 
time  steps.  Final  averages  give  an  indication  of  the  state  of 
the  copper  over  the  entire  undisturbed  period. 
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Alternative  Configurations 

The  Baskes  and  Foiles  code  allowed  further  operator  input 
in  order  to  specify  the  type  of  problem  to  be  solved. 

Different  parameters  such  as  volume  or  total  energy  can  be 
held  constant  depending  on  the  atomic  ensemble  to  be  generated 
by  the  experiment.  One  choice  is  between  fixed  or  dynamic 
crystal  boundaries.  In  order  to  use  changes  in  volume  as  an 
indication  of  the  state  of  the  copper  crystal,  the  researcher 
selected  dynamic  boundaries  for  all  computer  runs. 

Another  choice  was  between  a  constant  energy  problem  and  a 
temperature  fixed  problem.  In  the  constant  energy  situation, 
the  internal  conditions  of  the  crystal  reflect  the  reactions  of 
the  crystal  as  if  no  energy  is  allowed  to  enter  or  leave  the 
simulation.  In  the  temperature  fixed  situation  an  external 
temperature  is  applied  to  the  crystal  and  energy  can  flow  into 
or  out  of  the  crystal  and  into  the  surrounding  environment. 

All  computer  runs  to  gather  data  on  the  state  of  copper  were 
done  with  the  temperature  fixed  so  that  the  crystal's  internal 
energy  could  be  affected  and  so  that  the  state  could  change. 


Time  Saving  Steps 


In  order  to  save  computer  time  and  still  obtain  useful 


answers,  time  saving  steps  are  available.  Rather  than  compute 
the  effect  of  every  atom  on  every  other  atom,  for  example,  the 
computer  can  develop  a  neighbor  list  for  each  atom.  Neighbors 
are  those  atoms  located  within  a  specified  cutoff  distance, 
also  called  a  neighbor  interaction  sphere  as  in  Figure  2-4. 


Atoms  outside  the  cutoff  are  assumed  to  exert  negligible  force 
on  the  individual  atom  under  evaluation;  hence,  by  definition 
within  the  computer,  atoms  outside  the  cutoff  radius  do  not 
directly  interact. 

Either  the  original  neighbor  list  can  be  used  through  the 
rest  of  the  computer  run  or  the  list  can  be  revised  when  any 
atom  is  displaced  more  than  a  threshold  distance.  Both  the 
cutoff  distance  and  the  threshold  distance  can  be  specified  by 
the  operator  or  can  be  calculated  by  the  computer. 

In  order  to  avoid  repeating  calculations,  steps  can  be 
eliminated  by  applying  Newton's  Third  Law  of  Motion:  Every 
action  has  an  equal  and  opposite  reaction  (11:66).  This  means 
that  every  atom  affects  every  other  atom  with  the  same  force 
in  the  opposite  direction  and  eliminates  the  need  to  calculate 
the  forces  between  atoms  j  and  i  if  the  forces  between  atoms  i 
and  j  are  already  known. 

An  advanced  time  saving  step  is  to  use  a  different 
computational  approach,  a  parallel  processing  computer.  This 
approach  is  discussed  in  the  next  chapter. 
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1 1 .1 .  PARALLEL  PROCESSING 


In  order  to  expedite  program  execution,  the  molecular 
dynamics  program  was  modified  to  run  on  a  parallel  processing 
computer.  This  chapter  discusses  the  definition  of  parallel 
processing,  how  certain  molecular  dynamics  calculations  are 
well  suited  for  solution  by  parallel  processing  techniques, 
and  how  the  specific  embedded  atom  program  was  modified  to  run 
in  parallel. 

Parallel  Processing  Computers 

The  conventional  approach  within  a  computer  is  serial 
processing,  the  practice  of  a  single  computer  working  through 
a  problem  one  step  at  a  time.  Parallel  processing  is  often 
mentioned  as  one  of  the  next  great  advances  in  data  processing 
because  may  allow  many  steps  to  be  done  at  once  (5:13).  The 
basic  idea  is  to  have  a  controlling  computer,  also  called  a 
host,  manage  several  computers,  called  nodes,  that  can 
communicate  with  the  controller  and  with  each  other.  The  node 
computers  can  be  given  the  necessary  information,  complete 
their  given  tasks,  and  send  the  results  to  other  nodes  or  to 
the  host.  With  this  arrangement,  each  node  can  work  part  of  a 
larger  problem  and  many  nodes  working  at  once  can  reduce  the 
time  needed  to  reach  a  solution. 


Modification  of  any  program  from  a  serial  to  a  parallel 
approach  requires  the  programmer  to  select  the  processes  that 
can  be  solved  in  parallel  and  the  communications  arrangements 


part  of  the  array  and  communicates  its  results.  Many 
host/node  communication  arrangements  are  possible,  ranging 


need  to  store  either  fewer  individual  neighbor  arrays  or  one 
much  shorter  neighbor  array.  For  example,  the  node  which 
works  atoms  eight  through  fifteen  stores  the  neighbor  lists 
for  only  those  atoms.  That  node  does  not  need  to  waste  time 
or  memory  by  calculating  and  listing  the  neighbors  of  atoms 
outside  its  assigned  area. 

Force  calculations  can  be  made  much  faster  by  the  use  of  a 
similar  approach.  Each  node  can  calculate  forces  on  an 
appropriate  share  of  the  atoms  in  the  lattice  and  report  its 
results  back  to  the  host.  The  host  combines  those  results  in 
order  to  continue  the  calculations  for  the  particular  time 
step. 


The  Embedded-Atom  Solution  in  Parallel 


An  embedded-atom  program  modification  which  uses  parallel 
processing  and  which  gives  complete  solutions  follows  the 
basic  approach  outlined  above.  The  communications  plan  is  to 
have  the  host  send  identical  information  to  all  nodes.  Each 
node  first  assigns  itself  several  atoms  within  the  copper 
crystal.  Each  individual  node  finds  the  neighbors  of  its 
assigned  atoms,  then  determines  the  local  electron  density  for 
each  assigned  atomic  position.  The  electron  density  is  used 
to  find  the  embedding  energy  of  each  atom,  and  finally  each 
node  calculates  the  core  repulsion  between  its  assigned  atoms 
and  all  other  atoms  on  its  neighbor  list.  The  results  are 
passed  back  to  the  host  and  the  nodes  are  ready  to  begin 
again. 
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The  obvious  advantage  of  the  parallel  approach  is  that  the 
work  is  divided  among  several  processors  and  the  forces  and 
neighbors  for  many  atoms  are  solved  at  once.  This  proccess  is 
illustrated  by  Figure  3-3.  Adding  nodes  does  not  divide  the 
total  computer  time  by  the  number  of  nodes,  however. 
Communications  between  host  and  nodes  becomes  time  consuming 
as  only  one  node  can  report  back  to  the  host  at  one  time  and 
as  the  host  computer  must  reinitialize  the  arrays  used  for 
communications  after  the  reports  are  received.  The  potential 
advantages  of  parallel  processing  become  a  trade  off  between 
time  saved  by  simultaneous  operations  and  time  lost  to 
increased  communications  reguirements . 

Unfortunately,  also,  the  parallel  processing  computer  was 
not  continuously  available  during  this  research  period  and  the 
program  modification  process  proved  to  be  quite  time 
consuming.  Because  the  iPSC  cannot  run  in  a  background  mode, 
parallel  computer  runs  were  restricted  to  provide  validation 
and  timing  information  rather  than  provide  scientific  data  on 
the  copper  crystal.  Timing  advantages  are  discussed  in 
Chapter  Four  --  Experiments  and  Results.  The  scientific  data 
used  in  the  investigation  of  the  equation  of  state  was 
generated  on  the  ELIXI  system  with  its  sequential  processor. 


study  investigates  the  phase  transitions  of  copper  as  external 
pressure  was  held  constant  and  a  selected  external  temperature 
was  applied  and  held  constant  through  each  entire  computer 
run.  The  resultant  changes  in  internal  pressure,  volume, 
internal  temperature,  and  potential  energy  were  observed.  A 
change  from  the  solid  to  a  liquid  state  was  observed  when  the 
volume,  linear  displacement,  and  potential  energy,  as  plotted 
against  the  independent  variable  of  internal  temperature,  made 
a  discontinuous  jump  from  what  was  otherwise  a  normal  pattern 
of  expansion  and  contraction  driven  by  temperature. 

A  risk  of  using  discontinuities  in  volume  and  energy  is 
that  the  transition  may  not  be  readily  apparent  if  the  crystal 
is  caught  in  a  partially  melted/frozen  state  where  the  average 
values  fall  between  the  values  for  liquid  or  solid.  Although 
those  points  are  hard  to  analyze  at  first,  such  ambiguous  data 
points  provide  opportunities  for  further  study  because  the 
final  state  of  the  atoms  of  that  run  can  be  used  to  begin 
another  run  for  a  longer  time  or  at  a  slightly  different 
temperature.  The  combined  information  from  such  multiple  runs 
can  provide  increasingly  specific  information  on  melt  and 
freeze  temperatures. 

Another  problem  of  the  molecular  dynamics  approach  is  that 
artificial  superheating  and  supercooling  generally  cause  one 
to  miss  the  actual  melt  or  freeze  points.  Several  aspects  of 
molecular  dynamics  combine  to  make  this  a  possibility  that 
must  be  recognized  and  not  allowed  to  invalidate  otherwise 
good  results  from  computer  studies.  Because  the  crystals 


simulated  in  the  computer  are  necessarily  small,  time  steps 
are  too  small  to  be  directly  observed,  and  the  mathematics 
allows  pure  crystals  in  perfectly  defined  conditions,  some 
results  not  readily  observable  in  the  laboratory  must  be 
accounted  for.  Superheating  was  a  problem  in  the  first 
attempt  to  melt  the  copper,  and  supercooling  occurred 
consistently  in  the  cooling/freezing  experiments. 

The  entire  series  of  computer  experiments  on  the  ELIXI 
serial  processing  computer  used  a  copper  crystal  of  256  atoms. 
The  initial  run  began  at  293°  Kelvin  (K),  room  temperature, 
and  with  a  lattice  constant  of  3.165A,  in  agreement  with 
laboratory  data  (19:12).  Next,  the  copper  crystal  was  heated 
in  order  to  find  its  melting  point.  After  the  copper  melted, 
the  crystal  was  cooled  until  it  froze.  The  frozen  crystal  was 
then  reheated  to  see  if  the  melt  point  changed,  and  finally 
the  remelted  copper  was  again  cooled  in  an  attempt  to  repeat 
the  freezing.  All  runs  applied  an  external  temperature  which 
varied  from  run  to  run  but  was  held  steady  throughout  each 
run,  and  a  constant  external  pressure  of  zero  bars.  The 
internal  temperature,  energy,  and  volume  of  the  crystal  was 
recorded  and  averaged  over  several  thousand  time  steps. 

Validation  --  Copper  at  Room  Temperature 

The  first  computer  run  used  the  structure  of  an  fee  crystal 


and  the  lattice  constant  of  copper  to  build  the  matricies  that 
describe  the  positions  and  velocities  of  every  atom  of  a 
copper  crystal.  Each  atom  was  assigned  random  motions  with  a 


Boltzmann  distribution  based  on  293°  Kelvin,  the  external 
temperature  was  held  at  293°  K,  and  the  system  was  given  a 
standard  0.30  ps  (300  time  steps)  to  allow  the  system  to 
equilibrate.  The  crystal  was  given  flexible  boundaries  and 
the  observable  parameters  were  averaged  over  two  thousand  time 
steps  (2.0  ps).  The  average  volume  was  used  to  compute  the 
average  lattice  constant  and  that  distance  was  compared  to  the 
lattice  constant  of  copper  under  laboratory  conditions. 

The  average  volume  found  by  the  computer  was  was 
3067.86  A3  (+  0.079  A3)  which  gives  an  average  lattice 
constant  of : 

( 3067 . 86  )  3 / 3 

-  =  3.633  A  (±  0.00003  A) 

4 

which  is  within  0.50%  of  the  accepted  lattice  constant  of 
3.615  A  at  room  temperature  (20:12). 

With  the  computer  results  for  conditions  at  room 
temperature  in  agreement  with  actual  laboratory  data,  the 
volume  of  3067.86  A3  at  293°  was  used  as  the  baseline  to 
compare  volumes  at  other  temperatures.  Three  dependent 
parameters  --  volume,  linear  expansion,  and  potential  energy 
--  are  used  in  the  following  sections.  In  order  to  show  the 
change  in  the  state  of  the  copper  crystal,  the  three  dependent 
parameters  are  plotted  against  average  internal  temperature 
for  each  computer  run.  Discontinuous  changes  in  the  three 
parameters  indicate  a  change  of  state  and  all  three  plots  must 
agree  for  a  change  of  state  to  be  unambiguous. 


Determination  of  the  Equation  of  State 
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The  First  Series  of  Heating  Runs 

The  first  task  was  to  determine  if  the  copper  would 
melt  at  the  accepted  temperature  of  1357.6°  K  (18:77).  This 
was  done  by  first  building  a  computer  model  of  a  crystal  at 
293°  K,  heating  it  to  300°  K  and  holding  it  there  for  2.3  ps 
while  using  the  internal  parameters  of  the  last  2.0  ps  to 
calculate  overall  average  values.  The  crystal  was  then  heated 
to  600°  K,  then  900°  K,  1200°  K,  1500°  K,  and  finally  1800°  K 
in  separate  computer  runs  before  it  melted.  This  melting 
point  is  higher  than  found  in  the  laboratory  and  it  originally 
put  the  entire  approach  in  doubt. 

A  Second  Form  of  Validation 

A  second  means  of  validation  was  then  applied  by 
comparing  the  percentage  of  linear  expansion  calculated  in  the 
simulation  against  accepted  tabulations  (18:77).  The  results 
are  plotted  on  Figure  4-1.  The  original  data  points  at  1200° 

K  and  1500°  K  were  higher  than  the  prediction,  but  after  the 
crystal  was  held  at  each  temperature  for  5.3  ps,  the  linear 
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distance  lowered  to  within  excellent  agreement  with  the 
expansion  predictions.  These  points  are  also  plotted  on 
Figure  4-1  and  they  help  confirm  that  the  computer  gives 
reliable  information. 
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Figure  4-1, 


Computer  Calculated  and  Expected  Linear 
Expansion  vs.  Temperature 

(Expected  values  from  18:77) 
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Melting  After  Superheating 

The  high  melting  temperature  of  the  first  series  of 
runs  is,  therefore,  an  apparent  example  of  superheating.  The 
results  are  shown  on  figures  4-2,  4-3,  and  4-4,  and  they  give 
unambiguous  agreement  that  the  change  of  state  occurred 
between  1500°  K  and  1800°  K.  All  three  graphs  include  data 
points  at  1950°  K  and  2100°  K  to  establish  trend  lines  for  the 
liquid  state. 

In  an  attempt  to  break  the  superheating  and  melt  the 
crystal,  the  copper  was  held  at  1500°  for  5.3  ps,  but  no  change 
of  state  was  observed.  Heating  the  crystal  to  1550°  and  1600° 
also  did  not  cause  a  change  of  state.  Temperature  changes  from 
1500°  to  1650°  and  from  1600°  to  1650°  K  melted  the  copper  as 
shown  on  Figures  4-5,  4-6,  and  4-7. 

To  summarize  the  first  series  of  heating  the  copper 
crystal,  the  simulations  resulted  in  superheating  the  crystal 
past  the  accepted  laboratory  melting  point  until  it  melted 
between  1600°  K  and  1650°  K. 

Freezing  the  Melted  Crystal 

Obviously  an  actual  liquid  would  not  hold  its  shape; 
a  melted  copper  crystal  would  form  a  molten  puddle  and  its 
atoms  would  loose  the  order  imposed  by  the  crystal  lattice. 

The  computer  simulation,  however,  allows  the  lattice  to 
maintain  its  order  even  though  the  crystal  has  expanded  and 
the  increased  distances  between  atoms  allows  the  interatomic 
forces  to  decline  enough  to  allow  the  change  from  solid  to 
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Figure  4-2.  Volume  of  the  Copper  Crystal  as  it  is  Melting 
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Figure  4-8.  Volume  of  the  Copper  Crystal  as  it  Freezes 
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Figure  4-10.  Linear  Expansion  of  the  Copper  Crystal  as  it 
Freezes 
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Melting  the  Crystal  Again 

The  refrozen  crystal  at  1200°  K  was  again  heated  in  an 
attempt  to  either  repeat  the  superheating  or  find  a  more 
realistic  melt  temperature.  The  first  run  of  this  series,  to 
1500°  K  melted  the  crystal.  The  two  temperatures  of  1200°  K 
and  1500°  K  bracket  the  known  melt  temperature  of  1357.6°  K. 

The  next  run  used  the  same  start  point  at  1200°  K  and  heated 
the  crystal  to  1400°  K.  The  crystal  melted  again.  Heating  the 
1200°  K  crystal  to  1350°  K  gave  a  data  point  slightly  above  the 
trend  line  for  the  crystalline  state,  and  holding  that  crystal 
at  1350°  K,  7.6°  K  below  the  known  melt  point  of  copper,  the 
internal  temperature  averaged  1347.78°  K  and  the  observable 
parameters  fell  between  the  trend  lines  for  solid  and  liquid. 

In  order  to  confirm  that  the  copper  was  transitioning  near 
the  accepted  laboratory  melt  point,  the  slightly  melted 
crystal  at  1350°  K  was  cooled  to  1300°  K,  where  it  remained 
solid.  These  results  bracket  the  melting  point  of  copper 
between  1300°  K  and  1400°  K  and  at  least  imply  that  the  melt 
point  is  very  near  1350°  K.  These  results  are  shown  in 
Figures  4-11,  4-12,  and  4-13. 

With  the  remelted  copper  at  1400°  K,  there  was  one  last 
attempt  to  freeze  it  near  the  known  freeze  temperature.  The 
crystal  was  cooled  to  1300°  K,  but  the  data  points  remained  along 
the  liquid  trend  lines.  The  crystal  had  supercooled  again. 

Based  on  these  results,  the  melting  temperature  of  copper 
predicted  by  the  simulation  is  1350°  K,  plus  or  minus  50°  K. 
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Figure  4-12. 


Potential  Energy  of  the  Copper  Crystal  at  the 
Actual  Melt  Point 

The  arrows  represent  how  the  crystal  reacts 
to  different  temperatures. 
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Diagnostic  devices  available  on  the  UNIX  operating  system 
show  that  from  96%  to  98%  of  the  central  processing  unit  (cpu) 


¥ 


time  needed  to  solve  a  molecular  dynamics  problem  on  a 
conventional  linear  processing  computer  is  spent  doing  the 
force  and  neighbor  calculations.  To  speed  the  solution 
process,  the  force  and  neighbor  subroutines  were  exported  to 
the  nodes  of  the  iPSC  parallel  processing  computer.  Test 
programs  were  run  using  different  numbers  of  nodes  in  order  to 
show  the  proportionate  time  advantages  of  dividing  the 
calculations  among  the  nodes. 

In  order  to  keep  the  timing  results  between  the  serial  and 
parallel  computers  comparable,  both  types  of  computers  worked 
with  identical  256  atom  copper  crystals  held  at  a  constant 
external  temperature  and  with  flexible  boundaries.  Programs 
which  make  comparable  numbers  and  types  of  force  and  neighbor 
calculations  were  run  on  both  types  of  computers. 

The  original  computer  code  used  special  timing  functions 
available  on  the  Cray  supercomputer.  Computers  based  on  the 
UNIX  and  EMBOS  operating  systems  are  not  able  to  measure  time 
with  the  same  sophistication.  Relatively  simple  functions 
were  used  to  measure  the  time  required  to  solve  the  entire 
atomistic  simulation  problem  and  to  measure  the  time  spent  in 
the  subroutines  exported  by  the  parallel  computing  process.  A 
necessary  assumption  for  this  section  is  that  the  timing 
routines  added  to  the  original  computer  code  gave  sufficiently 
accurate  timing  information. 
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A  simple  comparison  of  total  time  spent  on  each  program 
would  not  give  meaningful  data  because  the  ELIXI  and  iPSC 
computers  run  at  different  speeds,  and  because  total 
processing  time  depends  on  the  level  of  multitasking  which  is 
beyond  the  control  of  the  experimenter.  Therefore,  in  order 
to  establish  potential  time  savings,  the  critical  difference 
between  serial  and  parallel  computer  runs  was  the  proportion 
of  time  spent  in  the  subroutines  exported  to  the  parallel 
computer  nodes.  The  results  and  the  calculations  below 
establish  that  the  parallel  approach  has  the  potential  to  save 
considerable  computer  time. 

In  order  to  make  complete  and  accurate  comparisons  the 
serial  UNIX  programs  were  run  in  both  a  normal  and  an 
unfavorable  configuration  which  did  not  use  all  available  time 
saving  steps.  The  normal  configuration  updated  the  neighbor 
list  only  as  necessary.  Cumulative  time  spent  in  the  force 
and  neighbor  subroutines  was  96%  of  the  total  time  of  the 
computer  run.  An  unfavorable  configuration  was  run  where  the 
neighbor  list  was  updated  at  every  time  step.  The  force  and 
neighbor  subroutines  used  98%  of  the  total  computer  time. 

After  the  program  was  rewritten  for  the  parallel  computer, 
experiments  were  run  to  obtain  similar  timing  information. 
Unfortunately,  the  parallel  processing  modification  at  this 
time  defaults  into  a  mode  which  updates  the  neighbor  list  at 
every  time  step.  Corrections  to  this  problem  have  yet  to 
work.  Because  neighbor  list  recalculation  is  a  time  consuming 
step,  the  only  comparison  fully  substantiated  by  experiment  at 


V.  CONCLUSIONS AND  RECOMMENDATIONS 


Instances  of  superheating  and  supercooling  not  normally 
observed  in  the  laboratory  can  be  accepted  because  of  the 
small  size  of  the  crystal  and  because  of  the  short  time 
periods  --  two  to  five  picoseconds  --  which  can  be  observed  by 
the  computer  simulations.  When  the  crystal  was  given  enough 
time  to  react  to  outside  forces,  it  melted  or  stayed  solid  at 
the  accepted  laboratory  temperatures.  In  all,  despite  the 
limits  of  atomistic  simulation,  the  molecular  dynamics  program 
and  the  embedded-  atom  method  gave  results  that  encourage 
their  use  in  further  studies. 

Parallel  processing,  in  general,  still  looks  promising 
despite  the  researcher's  inability  to  use  all  of  the  available 
time  saving  steps.  Additional  work  on  the  program,  improved 
programming  techniques,  and  advanced  computer  hardware  can  all 
combine  to  make  the  needed  changes.  At  this  time,  computer 
reliability  and  availability,  problems  with  the  parallel 
operating  system,  and  the  limits  on  host/node  communications 
and  memory  sharing  all  degrade  the  utility  of  the  iPSC  as  a 
scientific  tool. 

A  central  problem  for  the  iPSC  at  the  Air  Force  Institute 
of  Technology  is  the  lack  of  node  multitasking;  only  one  user 
can  use  the  node  configuration  at  one  time.  This  lack, 
combined  with  the  multitasking  of  the  host  processor,  allows 
bottlenecks  and  delays,  and  means  that  the  atomistic 
simulation  experiments  must  be  kept  relatively  short  so  that 
other  operators  can  continue  their  research. 
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processing  time  and  increase  accuracy.  Communications 
requirements  between  the  host  and  the  nodes  can  be  reduced  by 
memory  sharing,  which  allows  all  the  processors  to  access  the 
same  memory  space  rather  than  send  all  information  through 
time  consuming  send  and  receive  subroutines  (7:72).  An 
operating  system  which  allows  multitasking  of  the  nodes  will 
greatly  improve  the  availability  of  the  parallel  processing 
computer.  With  that  improvement,  several  researchers  will  be 
able  to  run  their  programs  simultaneously,  as  now  possible  on 
the  serial  processing  computers. 

As  computer  capability  improves,  the  size  of  the  crystals 
simulated  by  the  computer  can  be  increased  and  this  advance 
should  allow  better  accuracy  and  should  avoid  the  problems 
observed  in  this  study  such  as  artificial  superheating  and 
supercooling.  These  possibilities  should  combine  with  the 
proven  ability  of  molecular  dynamics  and  the  embedded-atom 
method  to  make  atomistic  simulation  a  more  useful  tool  for  the 
study  of  matter. 
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APPENDIX  A 


CALCULATION  OF  THE  POTENTIAL  TIME  SAVINGS 


The  purpose  of  this  appendix  is  to  show  the  calculations 
of  the  time  saved  by  the  11%  time  advantage  discussed  in 
Chapters  4  and  5. 

Recall  that  the  atomistic  simulation  program  spent  96%  of 
its  time  doing  the  force  and  neighbor  calculations  in  a  serial 
computer  and  85%  of  its  time  doing  the  same  calculations  in  a 
parallel  processing  computer  with  32  nodes. 

A  typical  computer  experiment  of  2300  time  steps  ran  for  a 
total  of  1  hour  and  50  minutes  in  the  background  mode  in  the 
ELIXI  serial  processing  computer  (as  indicated  by  the  time 
at  the  beginning  and  end  of  a  print  out).  That  amounts  to  6600 
seconds  (sec)  of  total  wall  clock  time  spent  in  the  computer 
(But  much  less  than  6600  sec  of  CPU  time.  The  timing 
functions  were  not  used  to  give  total  CPU  time.  It  is  assumed 
that  multitasking  beyond  the  control  of  the  experimenter 
delayed  all  parts  of  the  program  equally. ) . 

Because  4%  (100%  -  96%  =4%)  of  the  computer  time  is  used 
to  do  all  calculations  beyond  the  force  and  neighbor 
subroutines,  the  total  time  spent  doing  the  rest  of  the 
program  is  264  sec  (6600  sec  x  4%  =  264  sec)  for  a  serial 
processor.  With  those  same  subroutines  exported  to  the  nodes 
of  a  parallel  processing  computer,  the  only  significant 
changes  to  the  computer  program  are  in  the  force  and  neighbor 
subroutines;  the  rest  of  the  program  remains  essentially 
unchanged  and  should,  therefore,  use  the  same  total  time  of 


264  seconds.  The  calculations  not  exported  then  take  up  15% 
(100%  -  85%  =  15%)  of  the  total  computer  time.  In  other 
words,  the  same  264  seconds  make  up  15%  of  the  computer  time 
used  by  a  parallel  processing  computer.  This  leads  to  the 
expectation  that  1760  sec  (264  sec/15%  =  1760  sec)  will  be  the 
total  wall  clock  time  required  for  a  parallel  processing 
computer  which  operates  at  the  same  processor  speed  as  the 
serial  computer. 

Thus,  the  11%  time  advantage  of  using  a  parallel 
processing  computer  to  do  the  force  and  neighbor  calculations 
allows  the  parallel  processor  to  do  the  entire  atomistic 
simulation  in  26.7%  (1760/6600  =  26.7%)  of  the  original  time 
for  a  time  savings  factor  of  3.75  (1/26.7%  =  3.75). 
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The  equation  of  state  of  copper  was  studied  by  using  computer 
simulations.  An  existing  molecular  dynamics  program  written  for  the 
Cray  was  modified  to  run  on  both  the  VAX  and  on  a  parallel  processing 
system.  The  program  applied  the  embedded-atom  method  to  determine 
interatomic  forces. 

Experiments  on  the  VAX  found  the  copper  melting  point  toQbe  very 
near  1350  K  which  is  close  to  the  known  melt  point  of  1357.6  K. 

Measurements  taken  on  both  types  of  computers  showed  that  parallel 
processing  offers  a  time  savings  factor  of  3.75  over  serial  processing 
systems . 
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